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The evaluation of the performance of adiabatic annealers is hindered by lack of efficient algorithms 
for simulating their behaviour. We exploit the analyticity of the standard model for the adiabatic 
quantum process to develop an efficient recursive method for its numerical simulation in case of both 
unitary and non-unitary evolution. Numerical simulations show distinctly different distributions for 
the most important figure of merit of adiabatic quantum computing — the success probability — 
in these two cases. 


I. INTRODUCTION 

The ongoing debate whether the devices D-Wave One and D-Wave Two, which contain hundreds of supercon¬ 
ducting flux qubits, demonstrate quantum or classical annealing m, is to a large extent due to the fact that their 
quantitative simulation is at the moment impossible. Characteristically, there was no controversy over quantum 
behaviour of an 8-qubit register , where a quantitative agreement was found between the experimental data and 
the theoretical predictions. Though one cannot expect to be able to model classically a large enough quantum 
system cni, the development of efficient classical algorithms could hopefully enable us to model quantum systems 
just big enough to take over. We can also exploit the idiosyncratic features of certain quantum devices, which 
simplify their modelling by classical means. Here we demonstrate such an algorithm for modelling linear adiabatic 
evolution, including the case when decoherence is present. This result depends on the analyticity of the final state 
of the system as a function of evolution time, which is rigorously proven. Our small-scale numerical realizations of 
this algorithm show that some features of quantum annealing are comparable to the behaviour of D-Wave machines. 
Since at this stage the simulations do not account for the specific physical conditions that influence the machine’s 
performance we stop short of drawing conclusions as to the strength of signatures of quantumness in the actual 
device. Nevertheless, we do hope that the tools developed here will enable a better understanding of this important 
problem. 

An adiabatic process evolves the state vector of a quantum annealer over time T along the trajectory s |'0(s)), 
where the parameter s G [0,1] is the reduced time s = t/T. The initial value |'0(O)) is the ground state of the initial 
Hamiltonian Hi. When the process stops at t = T the state |'0(1)) is considered as an approximant to the ground 
state of the custom design Hamiltonian Hj, which encodes the solution of an optimization problem of choice. We 
consider the standard linear case, when the evolution is driven by the time-dependent Hamiltonian 

H{s) = il-s)H, + sHf, (1) 


according to the Schrodinger equation 

^|V>(s))=-*(T/h)iT(s)|^(s)). (2) 

Let the quantum system be an Wqubit register. In order to fix notation we describe the eigenstates and the 
corresponding eigenvalues of H{s): 

H(s) \m; s) = Em{s) \m; s), where ifo(s) < i?i(s) < ... < i? 2 N_i(s). (3) 

However, this assumption is not necessary to derive the analytic results presented in next section. As we will see, 
our analysis remains valid in every case of bounded (in particular, finite) Hamiltonians Hi,Hf. 

The purpose of analysis and simulation of 0 is estimation of the success probability 

p = ||nv^(i)f= ^ IWi)|j;i)P, (4) 

where H denotes the projector onto the ground state space of Hf, e.g. for a non-degenerate ground state P = 
|(^/>(1)|0; l)p. (In fact, in most cases of interest to us the dimension of the ground state space d exceeds 1. This is 
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rigorously accounted for in all numerical results presented in this article.) P provides a measure of the accuracy 
with which the adiabatic process finds the solution of the underlying problem. For the Hamiltonians considered 
here, the instantaneous ground state |0;s) is non-degenerate for 0 < s < 1 (see Remarks), although the ground 
state oi Hf may be highly degenerate. In a practical optimisation problem any of these degenerate states would be 
an equally useful solution. While there is a unique vector in the ground state space adiabatically connected to the 
initial ground state, this will in general be a superposition of computational basis states (a random one of which 
will be seen in the readout) and not necessarily of interest in the context of optimisation. 

Having fixed T and the number of steps and the Hamiltonian, we repeatedly conduct the numerical experiment. 
Final Hamiltonians are drawn from an ensemble of binary optimization problems and the simulation outputs the 
probability of success P for each instance oi Hf. We record the outcome of the series of experiments in a histogram 
of P values. The histogram of the probability of success has long been adopted as a means for testing hypotheses 
about adiabatic computing, see e.g. mi and, more recently, [1]. 


Remarks. It is interesting to briefly consider the evolution of the ground state projector s i—Hg, where Hi = H 
is as above. Given that in the case at hand the dependence of H(s) on s is linear and a fortiori holomorphic, 
the theory of analytic perturbations applies and it gives us the following picture: A standard general result, see 
Theorem 1.8, p. 70, [12], tells us that all the eigenvalues Ek{s) as well as the corresponding eigenprojections, 
including Hg, will depend on s analytically in the entire complex plane. A priori these analytic functions may have 
algebraic type singularities at a set of isolated exceptional points in the complex plane, where the number of distinct 
eigenvalues is lower than in the generic case. More precisely, the eigenvalues as functions of s are holomorphic in 
simply-connected domains that exclude the exceptional points. At the exceptional points an eigenvalue function 
may have, although does not necessarily have to have, a branching point. Moreover, at a branching point the 
corresponding eigenprojection will have a pole. A priori, such a scenario, i.e. hitting a branching point where the 
eigenprojection has a singularity, might be construed as a phase transition of the first kind. From this point of view 
it is not clear a priori why a result such as Theorem |H.l should hold. (Note that the result can be reformulated 
for complex variable s stating that the solution of Q is an entire function of s.) 

If we restrict the variable s to the real line, the picture outlined above is considerably simplified. Indeed, in 
this case H{s) is a self-adjoint operator, and satisfies the assumptions of Theorem 6.1, p. 120, [ 12 ]. This result 
assures us that both the eigenvalues and the eigenprojections are holomorphic around every point of the real 
line. In other words, even if some Sg on the real line is an exceptional point, the projection will not have a pole 
there and none of the eigenvalue functions will branch. However, we emphasize that the dimension of the ground 
state may jump discontinuously at the exceptional point, e.g. consider Eq{s) with the corresponding projection 
Hg and Ei{s) with Hg. It may be that Eq{so) = Ei{so). Both these projections depend on s holomorphically 
but the projection onto the ground state of H{so) will be H),^ = Hg^ -|- Hg,,. This explains why the projections 
corresponding to the eigenvalues depend on s continuously, while the projection to the ground state may have 
isolated jump discontinuities. In particular, the latter scenario is bound to play out for at least one point sg S [0,1] 
when the ground state spaces of Hi and Hf have different dimensions. At this juncture it is interesting to invoke 
a recent result of Zhang et ah, II3|. which shows that in the particular case of Hamiltonians we work with in our 
numerical simulations, see (29)-(30), the ground state is non-degenerate for all s G [0,1). This implies that in all 
cases for which the ground state oi Hf is degenerate the minimal gap tends to zero as s approaches 1 from the left. 


II. CONSEQUENCES OF ANALYTICITY 


Since Hamiltonian 0 depends on parameter s linearly, the solution of Eq. 0 may be expected to depend on s 
analytically, at least for small s. Thus, we look for solutions in the form of a Taylor series 

1 p{s) = 1 p 0 + Sipi + S^V '2 + S^V'3 + ■ • ■ (5) 

In fact, we will demonstrate that the radius of convergence of series characterizing solutions of Eq. 0 is infinite, 
i.e. the series converge for all s, see Theorem H.l below. In order to simplify notation, let us define auxiliary 
operators 

A = -i{T/h)Hi and B = -i(T/h) {Hf - H,). (6) 


Thus, is equivalent to 


^ms)) = {A + Bs)\f^{s)). 

as 


( 7 ) 
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Below we will consider this equation for general bounded operators A and and only later draw conclusions about 
the specific case (|6| ) relevant to the application at hand. 

Substituting ([^ in Q we readily obtain 


IV'o) = IV'(O)) 

\ipi) = A\i;o) 

: ( 8 ) 
IV'n) = i (A|V'„-l) +^|V'n-2)) for n > 2 


We use this recurrence as the core of numerical schemas for simulation of the solutions of Q. It allows relatively 
efficient, as compared to ODE simulation, computation of consecutive terms and consecutive partial sums of the 
series ([^. Numerical experiment shows that this method, implemented directly as a MATLAB script, is typically 
faster than the standard ODE solver ‘ode45’ by a factor exceeding 3 for iV < 12 qubits. For N = 14, in which case 
the memory constraint plays a significant role, the speed up factor is approximately 60, reducing the computation 
time from just under one hour to under one minute. Moreover, as described in Section V, a careful implementation 
of an optimized code allowed us to simulate systems with up to fV = 24. 

A similar procedure can be applied if the interpolation is a polynomial of order m in s: in that case the three-term 
recurrence relation will be replaced by an m -I- 2-order recurrence. Indeed, this may also be possible for an analytic 
interpolating function of known Taylor series. However, to demonstrate the method we only discuss the simplest 
(linear) interpolation here. 

Note that evaluation of the success probability Q requires only the computation of 

p = \n {'ipo + iIji + ip2 + ■ ■ ■) P- 


In other words, the solution is found, immediately as it were, at the point s = 1 without the need of finding all 
the intermediate states |'i/'(s))- This is in stark contrast to the computation based on an application of an ODE 
solver to Eq. (§. The current method has been tested against our earlier study of the correlation between gap and 
success probability [13]; in that work we used a conventional (Dormand-Prince) solver. Note also that if series (§ 
is known to converge absolutely, then |V'(s)) automatically satisfies ([^. Therefore, the main issue at stake is the 
estimation of the radius of convergence 


Rc 


limsup W'tpT, 

V n—>-oo 


|l/n 


(9) 


Here || || denotes the £2 norm. In general Rc depends on the constituents of the process Hi,Hf and possibly even 
|■^/'(0)). From the point of view of simulations it is ideal to have Rc = 00 , which ensures that |'^(s)) given in ([^ is 
the solution of for all (reduced) times s = t/T. If on the other hand, Rc < 1, then the series cannot be used to 
estimate the success probability. But, as we shall see, the series ([^ converges everywhere. 

Let us introduce the following notation: 


a:=||A||, 6:=||H||. 


( 10 ) 


Here, || || is the operator norm. Throughout the article we assume 0 < a, 6 < 00 . For the particular case of 

S , we have a = (T/h) ||iLi||, b = (T/h) \\Hi — Hf\\. The Hamiltonians are nontrivial and bounded, e.g. finite 
nensional, and, moreover. Hi ^ Hf. We have the following result: 

Theorem II.1. Rc = 00 , i.e. series converges absolutely in the entire complex plane and its limit ip{s) for 

s G [0,oo) is the solution of & satisfying the initial condition ■0(0) = V'o- 

Proof. Recurrence Q readily implies 

\lpn) = ^Pnli’o), ( 11 ) 


where P„ = Pn(A,B) are operators defined via recurrence 

Po = I 
Pi = A 

■ ( 12 ) 
Pn+i = APn + nBPn-i for n > 1. 
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Note that we have ||Po|| = 1) Ill’ll! = cl, and ||P„_|_i|| = \\APn + nBPn-i\\ < a||P„|| +n6||P„_i|j. We wish to estimate 
the rate of growth of ||P„||. To this end let us consider an auxiliary scalar sequence (p„) defined via recurrence 


Po =1 
Pi = a 

Pn+I = apn + nbpn -1 for n > 1. 


(13) 


Since po = ||Po||,pi = ||Pi|| and grows at least as fast as ||Pn||, we clearly have 

\\Pu\\<Pn. (14) 

Next, we undertake to estimate the growth rate of (pn)- First, observe that Pn = p„(a, &) may be viewed as a 
polynomial in two variables, e.g. P 2 = + b, ps = + Zab, etc. It is easily seen that the general form of the 

polynomials is 


Pn = co[n]a"’+ ci[n]a"' ^6 + C2[n]a"' ‘^b^ P C3[n\aJ^ + ... = ^ Cfc[n] i 


n—2k iJi 


(15) 




where the coefficients Cfe[n] remain to be found. Note that the last term of the polynomial is either Cn/ 2 [n] 5”'^^ 
(when n is even) or C[„/2][n] (when n is odd), with [x] denoting the integer part of x. Moreover, it is easily 

seen that (131 implies 

co[n] = 1 


ci[n] = (^) 


(16) 


Ck[n] = Cfc[n - 1] + (n - l)cfe_i[n - 2] 
Using this, and applying induction, one readily obtains an explicit formula 


Cfe[n] = {2k- 


where {2k — 1)!! = 1 • 3 • 5 • ... {2k — 1). In light of this (15) yields 


n.\ ^—L 


k=0 


1 

nn-2k O 

kl{n — 2k)l 2^ 


Next, we make the following observation 


fc!(n —2fc)!> — ! for fc = 0,1, 2,..., [n/2]. 

- O - 


Indeed, we either have k > [ f] or else n — 2k > n — 2 [^] > [^]. In either case (19) follows trivially. 
Next, we obtain from (|18| and (19): 


Pn — £ 


k=0 


_1_ ^n—2k ^ 

k\{n-2k)\^ 2^ 


- TFTT®"il + 2 ^ + (2^) +(2^) +■•■ +(2^)' ^ 




As is well known, ([^] —)■ oo as n —)■ oo. Thus, recalling definition (|ll[), we obtain 


ll^nir/" < Ib^nir/'ll^oir^" <{-^Pn] llV'oir/" = 

uaJ V 


nl 


ni 


llr. 




1/2 


which by ([^ implies Rc = 00 . □ 


1 /” ^ 0, 


(17) 


(18) 


(19) 


( 20 ) 


( 21 ) 

□ 
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Remark 1. Note that when b < 2a^ (valid in the adiabatic regime) estimate (20 1 may be replaced by a more 
efficient one. Indeed in such a case 


1 + 


b 


+ 





If in addition a < I, then terms Pn/n\ diminish very fast which ensures fast convergence of series ([^ when s = I, 
also in the numerical sense. However, this is not of interest in the current problem as it represents the anti- 
adiabatic limit, as the number of spin flips during the evolution is of order a. In the special case of interest 
a = {T/h)\\Hi\\,b = [T/h)\\Hi — Hf\\ and the condition a < I&6 < 2a^ is equivalent to 


^ implies \\Hi - Hf\\ < 2||R,||. 

This assumption on T ensures the most efficient computation of series (§. However, we do not claim that this 
restriction delineates the only regimes in which computation is effective. 


III. THE MASTER EQUATION 


Now we will extend the above results to the evolution described by the master equation (see e.g. my 

±p = -y[H,p]+ LpL^ - iLtLp - IpLti. 

Here, p is the density matrix, H is the quantum system Hamiltonian and the Lindblad operator L accounts 
for decoherence due to the uncontrollable interaction of the quantum system with the environment. In some 
applications the master equation involves a finite number of different Lindblad operators. While for our purposes 
it is sufficient to have just one Lindblad, the results presented here can easily be extended to the more general case. 
From now on we will only consider the evolution of finite-dimensional systems. Furthermore, as before we assume 
that H takes the special form ([^ with s = t/T. Substituting variables (|^ as before, we obtain an equivalent form 
of the master equation 


— p = Ca[p\ + sCb[p\ 


( 22 ) 


with the shorthand notation 


Ca[p] ■■= [A, p]+T (^LpL'^ - ]^L^Lp - ]^pL^ L}j , and 


Cb[p] ■■= [B,p\. 


(23) 


Note that the evolution equation (22) generalizes the von Neumann equation 


—p(s) = -i{T/h) [H{s),p{s)]. 


(24) 


In the context of adiabatic computing, it is vital to realize that there is an essential difference between the two 
types of evolution. Namely, at least if certain general conditions on Hf are satisfied, the adiabatic theorem applies 
to the case (241 and ensures that the probability of success increases with increasing T and will tend to one as T 


tends to infinity. In contrast, there is no known generalization of the adiabatic theorem, extending it to the case 
of ( [2^ . In fact, in the course of numerical experimentations we have observed that for a relatively large L the 
probability of success may decrease with increasing time T. Figure 3 illustrates this phenomenon. This observation 
agrees with the conclusions of Ref. m that there is a finite time window for the operation of a quantum annealer 
interacting with the thermal bath. 

In order to describe the analytic properties of ( [^ we introduce the Hilbert space structure on the linear space 
oi K X K matrices (e.g. K = 2^ for an N qubit system). For X G we define the Hilbert-Schmidt norm 

= Tr(AXl)^/^. This norm endows with the Hilbert space structure. We denote the corresponding 

operator norms of Ca,C,b '■ -a simply by |1/1 a||, ||'Cb||. As is well known, the Hilbert-Schmidt norm 
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has the submultiplicative property ||^i3||//s < ||A ||//5 
estimates: 


B\\hs which, when applied to (23), yields the following 


\\CA\\<2\\A\\HS + 2T\\LrHs, ||£B||<2|iB|U5. 

Next, we wish to consider solutions of ( |2^ . To this end we introduce the Taylor series Ansatz 

p(s) = po + spi + P 2 + + • ■ • 


(25) 


(26) 


As before, we substitute (26) into (22). In this way we obtain the following description of the Taylor series 
coefficients: 


Pn — I Bb')PO j 

n\ 

where Qn = Qni^A, £b) is a polynomial in variables Ca and Cb defined by recursion: 

Qo = I 

Qi = £a 

I^n+l — BaQti '^BbQu—I ffil" n ^ 1. 


(27) 


(28) 


We observe that the proof of Theorem can be repeated almost verbatim with obvious modifications such as: 
replacing A with Ca, B with Cb, and Pn with Qn, and redefining a := ||£a|| and b := |1£b| 1. This yields the 
following result: 


Theorem III.l. Series with coefficients defined via (STj ) and converges absolutely in the entire complex 
plane. Moreover, its limit p(s) for s € [0,oo) is the solution of (22) satisfying the initial condition p(0) = po- 

Remark 2. Consider the unitary map U{s) defined by the adiabatic Schrodinger equation: 

d 


ds 


U{s) = —i{T/h) H{s) U{s), so that t7(s)|^/>(0)) = |V'(s))- 


Note that recurrence (12) applied in the special case — i.e. A = —i{T/H) Hi, B = —i{T/h) {Hf — Hi) — provides 
a constructive description or simulation of U(s) via the formula 

oo 

n—0 

In particular this description of the semigroup s i—> [/(s) may be used to simulate the evolution of mixed states 
(p4). Indeed, 


—p(s) = -i{T/h) [i7(s),p(s)] 


p{s) = U{s)piO)U{s)*. 


IV. ALGORITHM ANALYSIS 


We will utilize elements of the proof of Theorem 
with the following specific problem constituents: 


k 

fee i,2,...,v. 

(29) 

~ ^Ising ~ 

— Jki c/CTf, k,l C 1,2,N. 

k<l 

(30) 


11.1 to analyze the efficiency of the recurrence algorithm 


Here, Jy S { — 1,1} denote random variables with uniform distribution P{Jij = 1) = P{Jij = —1) = 1/2, and N 
denotes the number of qubits. 
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First, in order to estimate the efficiency of a recurrence one needs to establish a stopping condition. It is 
desirable to stop recurrence when the approximation error is sufficiently small. However, little is known a priori 
about the solution that is sought in Q and so the calculus methods do not seem to provide an easy estimate 
for the dependence of error on the degree of approximation n. As a substitute we propose the following stopping 
condition: 


Compute for fc = 1, 2,... n, where n is the smallest index such that < e. 


(31) 


Here, e is regarded as a preset measure of accuracy. In light of the proof of Theorem |H.1| we know that n defined 
by (31) is bound to be finite. Moreover, since the solution dt(s) is an analytic function in the complex plane this 
is expected to provide a good estimate of the remainder error when s is bounded, e.g. s ~ 1. 

Next, consider the particular case T = 1 (recall also h = 1). Let us assume that the Taylor series is centred at 
So- In such a case a = \\Hi + so{Hi — Hf)\\ = 0{\\Hf\\) = b = \\Hi — Hf\\ = 0{\\Hf\\) = 0{N‘^) which 

readily implies 


^ + - 


1 . 


1/2 


a^ + -b =0{N^). 


(32) 


It now follows from (211, (31) and the Stirling formula that 

1/2 


O 1 + 


-V 

2a?) 


(SO 


1/n 


e n 


1/3 


i.e. n-^0(Ai®) 


Now, let us consider T » 1. While Theorems m and ensure that the Taylor series defining solutions of 
Eqns. ([^ and (22) converge in the entire complex plane, due to the roundoff errors, the Taylor series need not 
converge for relatively large values of T as intermediate sums become exponentially large in the adiabatic regime. 
The largest term /n\ in the Taylor expansion of occurs when n « |a:| and is of the order of To obtain 
solutions in such cases we use the following observation. Defining r = s — sq we observe that V’('r) satisfies (1^ 
where s has been replaced by t and A replaced by A-\- SqB with the initial condition '!/’|t=o = ■*/'(so)- Similarly, pyr) 
satisfies (22) where s has been replaced by r and A replaced by A + SqB with the initial condition p|t=o = p(sq)- 
This enables one to conduct computation in stages by partitioning the s-interval [0,1] into smaller segments or 
subintervals. Practice indicates that the best course of action is to carry out the computation in ^ T subintervals. 
Since computation within each segment is based on a faster converging Taylor series we obtain convergence of the 
solution in the entire s-interval for larger values of T. 

Note that the estimate (32) holds for every subinterval, which ensures that every segment requires the compu¬ 
tation of 0{N^) coefficients. Finally, it is easily seen that the number of arithmetical operations necessary to find 
n Taylor coefficients is of rank 0(n2^). Thus, the overall number of arithmetical operations required to carry 
out the computation is of rank 0{T ■ n • 2^) = 0{T ■ ■ 2^) = 0{2^). This is a polynomial time algorithm 

when measured against the classical problem size 2^. Naturally, and not very surprisingly, it is exponential in the 
number of qubits. The important feature of this algorithm is that it scales linearly with time T. 

As an example, we consider the Landau-Zener Hamiltonian H{s) = (1 — 2s)ct^ -|- It is easily seen that its 


minimal gap is given by ifi(s) — Eq(s) = 2-\/A^ -|- (1 — 2s)2, with the minimum occurring at s = .5. Consider the 
solution 'I'(s) of ([^ with d>(0) prescribed as the ground state of H{0). For the sake of an example let us attempt to 
compute 'k(l) for^ = 1, T/h = 20. In this instance the Taylor polynomial centered at s = 0 will not give meaningful 
results: experiment shows that the algorithm outputs d>(l) = [—0.0455 — 0.1663/,—0.4112 -|- 1.1521/] • 10^^ after 
100 iterations. Moreover one does not obtain correct convergence by simply increasing the number of iterations. 
Indeed, the norm of the simulated 4/(1) will converge to the value of ^ 1.05 instead of the correct ||4'(1)|| = 1. 
This is due to rounding error in the large intermediate terms. Indeed, all computations are done in floating point 
arithmetic with precision in the order of 10“^®. This means that the difference between any two consecutive real 
numbers in the order of 10^^ is as large as 10“®, which seriously impedes computational precision. In such a 
case a split of the interval into two parts, 0 < s < 0.5 and 0.5 < s < 1, results in smaller intermediate values 
and ensures convergence to the correct limit. This may be verified by a nearly exact calculation obtained via an 
application of, say, 200, instalments, see Fig. 11. The figure also shows the expected nonmonotonic convergence to 
unity with increasing T, with oscillations determined by the time-averaged Larmor frequency. In addition, we note 
that a symbolic implementation of our method may be used to reproduce the Taylor series of the special function 
solution. However, numerical evaluation of special functions is frequently quite problematic, and so the calculation 
via instalments offers a more practicable approach even if the special function solutions were known in closed-form. 

















Finally, we have used this specific example to benchmark the accuracy of our numerical method (termed TPT 
for Target-Point-Taylor) by comparing it with the standard fourth-order Runge-Kutta method (RK or ‘ode45’), 
m Based on a simple test we have observed the following: The TPT method gives the value of 

4'(1) ~ [0.509629891598850 -b 0.766898007985489*, -0.226356412675608 - 0.317659555887512 i], 

and 

P ~ 0.999801214304354 

with just two steps computing 350 Taylor coefficients at each step. Re-computation based on 200 intermediate 
steps shows that the value of 4>(1) is accurate to 10 significant digits while the value of P is accurate to 9 significant 
digits. Secondly, we have attempted to reach maximal accuracy via the RK method. In order to obtain comparable 
accuracy (7 significant digits for 4>(1) and 9 significant digits for P), namely 

4'(1) ~ [0.509629851149443 -b 0.766898035396728*, -0.226356395949650 - 0.317659566414614 i], 

and 

P ~ 0.999801214234416 

we were forced to use a very fine time-discretization ds = 1.86 x 10“® as larger time-steps lead to greater error. 
This example suggests that for comparable accuracy the ratio of the number of steps required by RK to those 
required by TPT exceeds 5.4 x 10® : 700. As mentioned above, the ODE-solver based methods are impracticable 
for larger systems due to long computation time. 

Remark 3. While the formal analysis above suggests that the constant of proportionality between the number of 
Taylor coefficients n and may be very large numerical evidence suggests that a rigid a priori number, 

say n = 50...90 is typically sufficient to ensure numerical convergence. However, for a small fraction of a percentage 
of Ising Hamiltonians there is no convergence under these conditions. Nevertheless, in the type of experiment 
considered in this article those cases are statistically insignificant and may be discarded. 

Remark 4. We note that the recurrence algorithm does not require that all the Taylor series coefficient terms 
be kept in the memory. Indeed, in the case of unitary evolution with the dynamic variable being a vector 4* of 
length 2^, the memory requirement is of rank 3 • 2^, regardless of the number of coefhcients n that are kept in the 
approximation. Indeed, we only need to keep in the memory the variable 4*, which is continually updated, and the 
two last Taylor coefficients, and as k runs from 2 to n. For a similar reason the memory requirement in 

the non-unitary case, when the dynamic variable is p of size 2^-^, is of rank 3 • 2^^. 


V. NUMERICAL EXPERIMENTS 


We have conducted numerical experiments for both a pure adiabatic process (§ and the process with dissipation 
(22) using matrix Hamiltonian as in (29) and ([30|). The choice of values for J^’s is generated via MATLAB’s 


generic random number generators. It is done independently for every run of the experiments discussed below. 
Note that \\Hi\\ = N and ||i7/|| < N{N — I)/2. 

Following the discussion on Section IV, we divide the s-interval in ~ T subintervals. Even though convergence 
in each subinterval is faster, the computational time is generally increased by a constant factor not exceeding the 
number of stages. 

The numerical simulations were implemented using compiled MATLAB code. The running time was lessened by 
using sparse matrices and refactoring the formulas ([^ and (|^ to reduce the number of matrix-vector calculations 
in favour of vector-vector calculations. The construction of Hi and Hf was made much more efficient by exploiting 
the structure and symmetries of the matrices. Memory requirements were also lessened by taking advantage of the 
sparseness and symmetry in the matrix computations. As a result, N up to 24 becomes computable on a modest 
distributed-computing server. 

Appendix H provides the source code of our current MATLAB implementation of the algorithm. The Main 
function executes the simulation the indicated number of times, using a sparse matrix H_i (the Heisenberg matrix 
representation of the initial Hamiltonian) and a different sparse vector H_f (the diagonal of the Random Ising Hamil¬ 
tonian in the Heisenberg matrix representation) at each time. Functions Initialfficuniltonian and Raindom_Ising 
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create those sparse structures using fractal-like iterations. Finally, function Taylor_Installments calculates the 
Taylor series. 

The histograms presented in Appendix I were produced using up to 180 cores of 2.0GHz Intel E5-2640L Xeon 
processors distributed on the compute cluster PLATO of the University of Saskatchewan. To illustrate the per¬ 
formance of the implemented algorithm, a total of 21780 individual simulations were ran for A^ = 20 in about 14 
hours, and using less than 5GB of memory per simulation. 

Figures 9 and 10 below illustrate the evolution of the compute time per simulation as a function of the number 
of qubits. An exponential increase of the compute time for large numbers of qubits can be observed. We expect 
this to be the case for any number of qubits and, in consequence, compute times for simulations of larger systems 
can be estimated. The fact that this exponential tendency is not verified when the number of qubits is too small, is 
certainly due to the fact that MATLAB, and the server itself, requires some time for internal processes (other than 
the execution of the algorithm). It appears that this time is much larger than the execution time of the algorithm 
for these cases when the number of qubits is too small. 

The crucial limiting factor in simulation of larger systems is memory capacity. However, early stage evidence 
shows that further algorithm optimizations, specialized to STB RAM server architectures, enables simulations of 
size up to 36 qubits. A detailed account of this ongoing work goes outside the scope of this article and will be 
reported elsewhere. 


A. Unitary adiabatic evolution: Schrodinger’s equation 

We use a schema based on recurrence (|^. Furthermore, we set |'(/'o) to be the ground state of Hi. We note that 
Hf is exponentially degenerate: the Hilbert space is 2^-dimensional but there are only 0{N) distinct eigenvalues, 
which are separated by multiples of 2. In particular, the Hamiltonian H{s) is invariant under a global spin flip for 
all time, so adiabatic evolution occurs in the subspace goanned by equal superpositions of computational states and 
the corresponding spin-flipped states, {(|a6c- • ■) -I- \abc- ■ •)) Let H denote the projector onto the ground 

state space of Hf. (Note that H is computed numerically, which presents no difficulties.) The probability of success 
is defined as P = ||n'!/)(l)|p. 

The results from simulation are presented via probability of success histograms. Fig. 1, and Figs. 4-8. Of note 
is the ragged landscape feature of the histograms. This can be understood qualitatively in the following way. The 
success probability is largely determined by the low-lying eigenvalue structure. The aforementioned degeneracy 
reduces the number of qualitatively different classes and we conjecture that the peaks correspond to different 
classes. We conjecture a connection between the low-P tail of the distribution of probability of success and the 
existence of optimization problems that are inherently hard. Since the optimization problem at hand is known to 
be NP hard, m, we expect this phenomenon to persist for larger systems. 


B. Non-unitary adiabatic evolution: Master equation with one Lindblad operator 

In this setting the probability of success is defined as P = Tr (np(l)) where, again, H is the projector onto the 
ground state space of Hf. In order to compute p(l) we use a schema based on recurrence (in notation of Subsection 

m]): 

Po = IV'o)(V'o|! where ipo denotes the ground state of Hi 
Pi = 

: (33) 

Pn+I = iPn + CbPti-i) for n > I. 

For the particular experiments presented in Figure 2 we took L = 0.1a where a is the operator, which in the 
eigenbasis of Hf assumes the off-diagonal form 


' 0 1 0 0 0 ... 
0 0 0 0 ... 
0 0 0 0 ... 


( 34 ) 
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The basis states are in order of increasing energy, with basis states in the degenerate subspaces listed in a particular 
order. The effect of this i7y-dependent Lindblad operator, in isolation, would be to relax the system to the selected 
ground state, thereby breaking the degeneracy. However, since it has a highly non-local action on the qubits, its 
interaction with H{s) disturbs the adiabatic evolution and destroys the fine structure of the probability histogram, 
see Fig. 3. 
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Appendix I: Figures 


10 ‘' Histogram: success probability for Schrodinger 8 qubits. ~500,000 runs 



Figure 1. Results obtained via computation of 
Taylor series solutions of the 8 qubit adiabatic 
Schrodinger equation with T = 4. The histogram 
displays probability of success outcomes for 500, 000 
random Ising Hamiltonians Hf sorted into 32 bins. 


10 “ Histogram: success probability for 8 qubits, T=4, L=.1a. ~100,000 runs 
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Figure 2. Results obtained via computation of Tay¬ 
lor series solutions of the 8 qubit master equation 
with the Lindblad operator L = 0.1a and T = 4. 
The probabilities of success outcomes obtained for 
100,000 random Ising Hamiltonians Flf are sorted 
into 32 bins. 


(a) Adiabatic Schrodinger evolution 
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(b) Adiabatic master equation evolution with L = .1* a 
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(c) Adiabatic master equation evolution with L = .05* a 

9 

10 

- 

0 

2 3 4 5 6 7 8 

(d) Adiabatic master equation evolution with L = .02*a 
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Figure 3. Dependence of the probability of success 
on T for solutions of the 8 qubit Schrodinger equa¬ 
tion (a) and the master equations with the Lind¬ 
blad operators L = 0.1a (b), L = 0.05a (c), and 
L — 0.02a (d). Note the inapplicability of the adia¬ 
batic theorem in (b) and (c), as well as convergence 
to the Schrodinger solution for diminishing norm of 
L. 


Success probability for Schrodinger 8 qubits, T = 10, -150,000 runs 



Figure 4- Schrodinger evolution as in Fig. 1 with 
8 qubits but with T = 10. The histogram may be 
described as multi-modal. 
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Histogram: success probability for pure Schrodinger 14 qubits 
T=10,19980 samples, with 128 bins 

6001- ^^^^^^^^-, 



Figure 5. Schrodinger evolution with 14 qubits and 
T = 10. The computation of these results took 
about 12 minutes using 90 cores. 


Histogram: success probability for pure Schrodinger 18 qubits 
T=10, 19980 samples, with 128 bins 

6001- ^^^^^^^^- 



Figure 7. Schrodinger evolution with 18 qubits and 
T = 10. Computation time was approximately 5 
hours and 20 minutes with 90 cores. 


Histogram: success probability for pure Schrodinger 16 qubits 
T=10, 19980 samples, with 128 bins 
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Figure 6. Schrodinger evolution with 16 qubits and 
T = 10. Computation time with 90 cores was 1 hour 
and 15 minutes approximately. 


Histogram: success probability for pure Schrodinger 20 qubits 
T=10, 19980 samples, with 128 bins 

6001- ^^^^^^^^- 



Figure 8. Schrodinger evolution with 20 qubits and 
T = 10. Using 180 cores, these results were com¬ 
puted in 14 hours. 
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Figure 9. Execution time (wall time) of the algo¬ 
rithm as a function of the number of qubits. Our 
current implementation running on MATLAB im¬ 
plies a system time larger than the time required 
by the algorithm when the number of qubits is too 
small. 


Execution time per simulation in a semi-log scale 



Figure 10. Execution time of the algorithm as a 
function of the number of qubits. A clear expo¬ 
nential behaviour is observed for larger number of 
qubits. The exponential function that fits the last 
three points (x=20, 22 and 24) is given by: y = ce™® 
where c ~ 6/rs, and m = 0.86. 
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91'P(s) 



Figure 11. The top two figures show a numerical solution of ([^ with 
T/h, = 20 and the Landau-Zener Hamiltonian H{s) = (1 — 2s)crz + Aca; 
obtained via condensing the Taylor series to 200 small steps. The star 
marked points are obtained via computation in just two instalments. The 
first (resp. second) figure shows the real (resp. imaginary) parts of the two 
vector components. The bottom graph shows oscillations of the probability 
of success P{T) as T/h varies between 20 and 50. 
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Appendix II: Source codes 


function P=Main(runs,N,T,nint,epsilon) 
dim=2" N; 

psi_0 = (1/sqrt(dim))*ones(dim/2,1); 
H_i=Initial_Hamiltonian(N-l); 
step=l/nint; 
for n=l:runs 

H_f=Random_Ising(N); 
ind=(H_f == min(H_f)); 
psi_in=psi_0; 
for k = l:nint 

psi_in = Taylor_Installements(H_i,H_f,T,psi_in, step,tol,k); 
end 

P[n] = 2*norm(psi_in(ind))" 2; 

end 


Function Main: Function that is being executed on each of the assigned cores on the cluster. Each core will obtains 
n =runs samples. Observe how the Taylor series is calculated for nint segments in the s-interval. N is the number 
of qubits, T the final time and tol the tolerance passed to the function Taylor_lnstallements below. 


function H=Random_lsing(d) 
dim=2" d; 

vals=zeros(l,dim/2); 
for k=l:d-l 

temp=zeros(l,dim/2~ k) ; 
for r=k+l:d 

mac=[ones(l,dim/2~ r) -ones(l,dim/2~ r)] ; 
for p=l:r-(k+l) 
mac=[mac mac]; 
end 

signo=2*randi(2,l)-3; 
temp=temp+signo*mac; 
end 

for p=l:k-l 

temp=[temp fliplr(temp)]; 
end 

vals=vals+temp; 

end 

H=sparse(vals’); 


Function Random_lsing: Using a fractal-like iteration, creates a sparse vector H of dimension 2“^ that represents 
the diagonal of an Ising Hamiltonian. Notice that this vector is palindrome, thus, only the first 2“^”^ values are 
returned. 
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function H=Initial_HamiltoniaLn(d) 

dim=2~ d; 
rows=[1 2]; 
cols= [2 1]; 
for k=l:d-l 

rows= [rows rows+2"k]; 
cols=[cols cols+2"k]; 
p=l:2-k; 

rows= [rows p p+2"k]; 
cols=[cols p+2~k p] ; 

end 

vals=-ones(l,size(rows,2)); 

H = sparseCrows, cols, vals, dim, dim); 


Function InitialJlamiltonian: Using a fractal-like iteration, creates a sparse matrix H (the initial Hamiltonian) 
of dimension 2'*. Notice that for N qubits we use Initial_Hamiltonian(N-l) due to the symmetry in H_f. 


function psi=Taylor_Installements (H_i, H_f _diag, T,psi_in, step, tol, k) 

Prepares Auxiliary Operators A and B of Eq. ([^ for the recurrence 
c = -li * T; 
d = (k - 1) * step; 

i_by_psi_in = H_i * psi_in - f lipud(psi_in) ; 
f_by_psi_in = H_f_diag .* psi_in; 
psi_n_min_2 = psi_in; 

psi_njmin_l = c * ((1-d) * i_by_psi_in + d * f _by_psi_in) ; 
psi = psijn_min_l * step + psi_n_min_2; 
nrm_cor = 1; 
n = 1; 

i_by_min_2 = i_by_psi_in; 
f_by_min_2 = f_by_psi_in; 

Implements the recurrence as in Eq. ([^ 
while (nrm_cor > tol) 
n = n+1; 

i_by_min_l = H_i * psi_n_min_l - f lipud(psi_n_min_l) ; 
f_by_min_l = H_f_diag .* psijn_min_l; 

psijn = (c/n)*((l-d) * i_by_min_l + d * f_byjmin_l + f_byjmin_2 - i_byjmin_2) ; 

cor = psi_n*step"n; 

nrm_cor = norm(cor) ; 

psi = psi + cor; 

psijn_min_2 = psi_n_min_l; 

psi_n_min_l = psi_n; 

i_by_min_2 = i_by_min_l; 

f_by_min_2 = f_by_min_l; 

end 


Function Taylor_lnstallements: Favouring vector-vector calculations, this function computes the Taylor series 
for the k-th segment of the s-interval. Notice that the calculation is stopped when the contribution at the next 
step is smaller than the tolerance tol. 




